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Abstract 

A damage model for the simulation of delamination propagation under high-cycle 
fatigue loading is proposed. The basis for the formulation is a cohesive law that links 
fracture and damage mechanics to establish the evolution of the damage variable in 
terms of the crack growth rate dA/dN. The damage state is obtained as a function 
of the loading conditions as well as the experimentally-determined coefficients of 
the Paris Law crack propagation rates for the material. It is shown that by using 
the constitutive fatigue damage model in a structural analysis, experimental re- 
sults can be reproduced without the need of additional model-specific curve-fitting 
parameters. 

Key words: Delamination, Fatigue, Cohesive elements. 


1 Introduction and motivation 


High-cyclc fatigue is a common cause of failure in aerospace structures. In 
laminated composite materials, the fatigue process involves several damage 
mechanisms that result in the degradation of the structure. One of the most 
important fatigue damage mechanisms is interlaminar damage (delamination). 

There are two basic approaches for the analysis of delamination under fatigue 
loading: Fracture Mechanics and Damage Mechanics. Fracture Mechanics re- 
lates the fatigue crack growth rate with the amplitude of the energy release 



rate and mode-ratio. In most studies, fatigue crack growth rates are described 
with the Paris Law [1,2]: 


d A _ ( AG \ m 

ON ~ C \G C ) 

where A is the crack area, N is the number of cycles, and the parameters 
C and m depend on the mode-ratio and must be determined experimentally. 
The cyclic variation of the energy release rate, AG, depends on the loading 
conditions, and G c is the fracture toughness of the material. Alternatively, 
the crack growth rate may be expressed in terms of the stress intensity factor 
range, A K, or the J-integral range A J [3,4], 

Damage Mechanics Models describe the loss of a material’s ability to carry 
loads by using one or more damage variables. In this paper, a Damage Me- 
chanics formulation is adopted in which fracture is represented by a cohesive 
zone model (CZM) that uses a single damage variable. The material separa- 
tion under cyclic loading is described by a constitutive equation formulated in 
the context of the thermodynamics of irreversible processes. Within an irre- 
versible cohesive zone model, material separation is described as a relationship 
between crack surface traction and separation across the crack. 

Linder general cyclic loading, the total damage is the sum of the damage 
caused by static or quasi-static loads and the damage that results from the 
cyclic loads. There are several models that extend cohesive laws for mono- 
tonic loading into forms suitable for cyclic loading. Yang et al. [5] developed a 
cohesive law that describes separately the unloading and reloading processes 
and creates a hysteresis loop between unloading and reloading paths. Roe 
and Siegmund [6] describe fatigue crack growth by incorporating a damage 
evolution equation for cyclic loading. Nguyen et al. [7] developed a cohesive 
zone model for cyclic loads in which the irreversible material degradation is 
represented as a loss of stiffness in the cohesive zone during the unloading 
portion of the load cycle. Similarly, Goyal-Singhal et al. extended the capabil- 
ity of their cohesive-decohesive constitutive model [8] to account for fatigue 
damage accumulation during unloading [9]. In all of these references, the fa- 
tigue damage accumulation is accounted for in a cycle-by-cycle analysis. For 
high-cycle fatigue, where the number of cycles considered is larger than 10 5 , 
a cycle-by-cycle analysis would be computationally intractable. 

For high-cycle fatigue, the damage evolution that results from cyclic loads 
is usually formulated as a function of the number of cycles and strains (or 
displacement jumps) [10-12], In these references, a damage evolution law ex- 
pressed in terms of the number of cycles is established a priori by adjusting 
several parameters through a trial-and-error calibration of the analysis. 
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In this paper, an approach is proposed whereby the evolution of damage 
derives from a Fracture Mechanics description of the fatigue crack growth 
rate. The approach is formulated using the cohesive zone model concept. A 
constitutive damage model previously developed by the authors for static or 
quasi-static loads [14] is enhanced to incorporate a damage evolution law for 
high-cycle fatigue. 

In the present model for fatigue damage, the evolution of the damage variable 
associated with cyclic loading is derived from a Fracture Mechanics description 
of the fatigue crack growth rate. Therefore, the proposed model is based on 
linking Fracture Mechanics and Damage Mechanics. The model relates dam- 
age accumulation to the number of load cycles while taking into account the 
loading conditions (load ratio, R , energy release rate, G, and fracture mode 
mixity). When used in a structural analysis, the model can simulate the de- 
pendence of the crack growth rate on these parameters. In addition to the 
Paris Law crack growth regime, the model also exhibits a threshold value for 
no growth as well as quasi-static tearing. The new fatigue damage model is 
implemented as a user- written element in ABAQUS [37] based on the cohesive 
finite element previously developed by the authors [14]. 


2 Cohesive zone model approach 


The CZM approach [15-17] is one of the most commonly used tools to simulate 
interfacial fracture. The CZM approach assumes that a cohesive damage zone 
develops near the tip of a crack. 

As mentioned above, cohesive damage zone models relate tractions, r, to dis- 
placement jumps, A, at an interface where a crack may occur. Damage ini- 
tiation is related to the interfacial strength, t°. When the area under the 
traction-displacement jump relation is equal to the fracture toughness, G c , 
the traction is reduced to zero and new crack surfaces are formed. If a linear 
softening law is used, the new crack surfaces are completely formed when the 
displacement jump is equal to, or greater than, the final displacement jump, 
A / (see Figure 1): 


r = (1 — d )t° 


A < A f 


t = 0 A> A f 


( 2 ) 


In the cohesive damage model, the damage variable d describes the density 
of microcracks of a representative element surface. Then, the damage variable 
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Fig. 1. Linear softening law for a cohesive zone model approach. 

can be interpreted as the ratio of the damaged area, A d , with respect to 
the area A e associated with the local discretization [18]. In the context of 
finite elements, the area A e represents the area of the element (or that of an 
integration point). Using the linear softening law represented in Figure 1, this 
ratio is a function of the energy dissipated during the damage process, 5, and 
of the critical energy release rate, G c . Using Equation (2), the damage variable 
can be expressed as the ratio between the current and the final displacement 
jumps: 


d = ^ = “_ = A 

yl e G c A f 


( 3 ) 


2. 1 Numerical representation of the CZM 


Cohesive finite elements have been developed to capture the initiation and 
propagation of delamination cracks [8,13,14,19-28]. As in most other cohesive 
element formulations, the constitutive law used in this paper is a bilinear 
relation between the tractions and the displacement jumps [13,14,26,29]. The 
bilinear cohesive law is similar to the softening law of the CZM but with an 
initial linear elastic response before damage initiation, as shown in Figure 2. 
This linear elastic part is defined using a penalty stiffness parameter, K , that 
ensures a stiff connection between the surfaces of the material discontinuity. 
The interfacial strength and the penalty stiffness define an onset displacement 
jump, A°, related to the initiation of damage. The equivalence between the 
constitutive equations of the physical cohesive zone model and the numerical 
constitutive equations is shown in Figure 2. 

The next section describes the kinematics and constitutive relation of cohesive 
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Fig. 2. Traction- displacement laws describing the physical (left) and numerical 
(right) constitutive equations of the CZM. 


zone models for quasi-static loading that are also the foundation of the fatigue 
damage model proposed in this paper. 


3 Kinematics and constitutive model for quasi-static loading 


The displacement jump across the interface [rq] , is obtained from the dis- 
placements of the points located on the top and bottom sides of the interface, 
uf and u ~ , respectively: 


M = uf - u t (4) 

where uf are the displacements with respect to a fixed Cartesian coordinate 
system. A co-rotational formulation is used to express the components of the 
displacement jumps with respect to the deformed interface. The coordinates 
Xi of the deformed interface can be written as [30]: 

Xi = + * (uf + u ~ ) (5) 

where X t are the coordinates of the undeformed interface. 

The components of the displacement jump tensor in the local coordinate sys- 
tem on the deformed interface, A m , are expressed in terms of the displacement 
held in global coordinates: 


A 


m 




( 6 ) 
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where Q„ u is the rotation tensor, defined in [13,14], 

The constitutive operator of the interface, D Jt , relates the element tractions, 
Tj, to the displacement jumps, Ap 


Tj = DjiAi (7) 

The constitutive relations of cohesive zone models must compute accurately 
the energy dissipated in the process of fracture. Under single-mode loading, 
controlled energy dissipation is achieved by ensuring that the area under the 
traction-displacement jump relation equals the corresponding fracture tough- 
ness. Under mixed-mode loading, a criterion established in terms of an inter- 
action between components of the energy release rates associated with each 
fracture mode is used to predict crack propagation. The formulation of the 
damage model previously proposed by the authors [14] is summarized in Ta- 
ble 1, where and A° are the free energy per unit surface of the damaged and 
undamaged interface, respectively. The function <5^- is the Kronecker delta, 
and the variable d is a scalar damage variable. The parameter A is the norm of 
the displacement jump tensor (also called equivalent displacement jump norm ) 
and it is used to compare different stages of the displacement jump state such 
that it is possible to define such concepts as ‘loading’, ‘unloading’ and ‘reload- 
ing’. The equivalent displacement jump is a non-negative, continuous scalar 
function defined as: 


^ — \J (A3) 2 + ( A s h ear ) 2 ( 8 ) 

where (•) is the MacAuley bracket defined as (x) = \ (x + |x|). The displace- 
ment jump in Mode I, i.e., normal to midplane is A3. The displacement jump 
tangent to the midplane, A shear , is computed with the Euclidean norm of the 
displacement jump in Mode 11 and Mode III: 


A shear ~ \[(Axf + (A2) 2 


( 9 ) 


The evolution of damage is defined by a suitable monotonic scalar function, 
G(-), ranging from 0 to 1. A damage consistency parameter, /i, is used to de- 
fine loading- unloading conditions according to the Kuhn- Tucker relations [31]. 
The consistency conditions ensure that damage evolution cannot occur during 
unloading or neutral loading. The damage threshold for the current time, t, is 
r x . The damage threshold is equal to the onset displacement jump, A° when 
there is no damage at the interface. When the interface is completely dam- 
aged, the damage threshold is equal to the final displacement jump, A A In the 
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Table 1 

Definition of the constitutive model. 


Free Energy 

if (A. d) = (1 - d) U° (At) - d U° (S 3i (— A 3 )) 

Constitutive equation 

n = = (f " d) SijKAj - d SijKStj (-A 3 ) 

Displacement jump norm 

A = \J (A 3) 2 + (Aghear) 2 

Damage criterion 

F (A 1 , r l ) := G (A 1 ) — G (U) < 0 Vt > 0 


r ( \ \ _ A^(A-A°) 
- X(Af-Ao) 

Evolution law 

j _ ;9F( \,r) _ ■ 0G(A) 
0 - d\ - n g\ 

Load/unload conditions 

fi > 0 ; F( A t ,r t )<0 ; /iF(A t ,r t ) = 0 


r l = max {A°, max s A s } 0 < s < t 


case of crack closure during load reversal, the constitutive model is designed 
to prevent interpenetration of the faces of the crack by restoring the normal 
penalty stiffness of the element even in the presence of damage. Further details 
regarding the damage model used here can be found in references [13,14], 

Under loading conditions, the damage variable can be obtained using the 
damage criterion in Table 1 as: 


, _ A'(A - A°) 
“ A(A/ - A°) 


( 10 ) 


In the numerical model of the CZM, the damage variable d represents a loss of 
stiffness and, therefore, it is not equivalent to d, the ratio between the damaged 
area, A d , with respect to the area of the element, A e , in Equation (3). Since d 
is equal to the ratio of the energy dissipated over the fracture toughness the 
damaged area ratio is related to the damage variable, d, as: 
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( 11 ) 


d = 


A* 

A e 


a 



(i-d) 


By solving Equation (10) for A and substituting into Equation (11), the dam- 
aged area ratio becomes: 


A d _ dA° 

A e ~ A/(l - d) + dA° 


( 12 ) 


4 Constitutive model for high-cycle fatigue 


The damage evolution that results from a general loading history can be con- 
sidered as the sum of the damage created by the quasi-static overloads and 
the damage created by the cyclic loads: 


d 


d 


static 


+ d 


cyclic 


(13) 


The Erst term in the right hand side of Equation (13) is obtained from the 
equations presented in previous section, while the second term has to be de- 
fined to account for cyclic loading. Using a Damage Mechanics framework, 
several authors have formulated the damage evolution that results from cyclic 
loads in terms of the number of cycles and of the strains (or displacement 
jumps) [10-12], These damage laws are model-specific and they are a func- 
tion of several parameters that have to be adjusted to calibrate the numerical 
model with experimental results, usually by trial and error. In contrast, the 
fatigue damage model formulated here is based on a Fracture Mechanics crack 
growth rate characterization which is achieved by linking Fracture Mechanics 
and Damage Mechanics: the evolution of the damage variable, d, is related 
with the crack growth rate, as follows: 


<9d <9d cL4 d 

dN = dA~ d ~dN 


(14) 


where A d is the damaged area, and ^ is the growth rate of the damaged area. 
The term Mf- is a material property that must be characterized experimentally 



for different loading conditions. The term 7^ can be obtained from Equation 

( 12 ): 


dd _ 1 [A f (1 — d) + dA°] 2 
d A d ~ A e AfA° 


(15) 


4-1 Determination of the growth rate of the damaged area as a function of 
the number of cycles 


In a degradation process involving cyclic loading, the damaged area grows as 
the number of cycles increase: after AN cycles, the damaged area ahead of 
the crack tip increases by AA^ as schematically represented in Figure 3. It 
can be assumed that the increase in the crack area A A is equivalent to the 
increase in the amount of damaged area. 




Crack tip 
position 



Fig. 3. Schematic representation of the equivalence between the increase in the dam- 
aged area and the crack growth. 


The increase in the damaged area along a crack front is equal to the increase 
in the damaged area of all of the elements ahead of the crack tip. Therefore, 
the crack growth rate can be assumed to be equal to the sum of the damaged 
area growth rates of all damaged elements ahead of the crack tip, that is, all 
elements in the cohesive zone: 
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( 16 ) 


dA 

dN 


E 

eeAcz 


dN 


where A d is the damaged area of one element and the term Acz is the area of 
the cohesive zone. Assuming that is the mean value of the damaged area 

dA e 

growth rate of the elements over the cohesive zone and assuming that the 
mean area of the elements in the cohesive zone is A e , the previous equation 
can be written as: 


d A _ dA d _ A C z dA d 

dN ~ ~dN ~ ~A^1)N 

eEAcz 


(17) 


where the ratio represents the number of elements in which the cohesive 
zone has been divided. In the context of finite elements, this ratio represents 
the number of elements that span the cohesive zone. Rearranging terms in 
Equation (17), the surface damage growth rate can be written as: 


dA d A e dA 
dN = AczM 


(18) 


4-2 Evolution of the damage variable under cyclic loading 

By introducing Equations (15) and (18) into Equation (14) the evolution of 
the damage variable as a function of the number of cycles can be written as: 


3d _ 1 (A7(l - d) + dA°) 2 dA 
M ~ A^ A/A° dN 


The area of the cohesive zone for pure Mode I can be estimated using Rice’s 
closed-form equation [32,33]: 


Tlcz — 


9t t E 3 G max 
b 32 (r°) 2 


( 20 ) 


where b is the width of the delamination front, and G max is taken as the 
maximum energy release rate in the loading cycle. E 3 is the Young’s modulus 
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of the bulk material in the direction perpendicular to the crack plane, and r° 
is the interfacial strength. 

4-3 Crack growth rate 

The crack growth rate under fatigue loading, |^, is a load and material- 
dependent characteristic that has been widely studied. The growth rate defined 
by the Paris Law given in Equation (1) represents crack propagation in region 
II of the typical pattern of the crack growth rate (see Figure 4). 



Normalized Energy Release Rate, 
log(G m 7Gc) 

Fig. 4. Typical crack growth rate regions. 


In region I, crack growth is not observed if the maximum energy release rate 
is smaller than the fatigue threshold of the energy release rate, G th- 
in. region III, the crack growth rate increases because the maximum energy 
release rate approaches the fracture toughness. Tearing fracture controls the 
crack growth rate in region III instead of fatigue propagation. 

The crack growth rate used in the fatigue damage model, Equation (19), 
is defined as a piecewise function defined as: 


dA 

dN 


C (jf) rn , G th < G max < G c 
0 , otherwise 


( 21 ) 
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where C, m and G c are material constants that depend on the mode-ratio. 
The maximum energy release rate G max and cyclic variation in the energy 
release rate AG used in the Paris Law rate equation can be computed using 
the constitutive law of the cohesive zone model (see Figure 5): 



A° l~ A 1 \ 


Fig. 5. Variation of the energy release rate. 
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It is clear from Equation (26) that the model accounts for variations in the 
load ratio. The higher the load ratio, the smaller the variation in the energy 
release, as shown in Figure 6. 


G 




Fig. 6. The load ratio effect is captured by the constitutive equations. The higher 
load ratios (R\ > R 2 ) the smaller AG (AG 2 > AG\). 


4-3.1 Mixed-mode loading 

The material parameters, C, m, Gth used in the crack growth rate expression 
(21) depend on the mode ratio. In Mode I, the crack growth rate parameters 
are Ci, mi, and Guh, and in Mode II, the crack growth rate parameters are Cn, 
mu , and Gn t h- Under mixed-mode, the crack growth rate parameters C , m, 
and G t h must be determined. I 11 this paper, the dependence of the parameters 
C and m on the mode ratio is assumed to be of the form [34]: 

logC' = logC' I +(^)logC7 m +(^ log^ ( 27 ) 


and 


m = 


mi + m m (§2) + (mu - m, - m,„) (ffj 


(28) 


where C m and m m are mode-ratio material parameters that must be deter- 
mined by curve-fitting experimental data. 


13 


The dependence of the energy release rate threshold is assumed to follow an 
expression similar to that introduced by Benzeggagh and Kenane [35] for the 
dependence of the Fracture Toughness with the mode-ratio: 


G th = G 


I th 


(G 


nth 



(29) 


where 772 is a material parameter obtained from a curve-fit of experimental 
results. 


4-4 Cycle jump strategy 


In a degradation process involving high-cycle fatigue, a cyclc-by-cyclc analy- 
sis becomes computationally intractable. Therefore, a cycle jump strategy is 
implemented in the finite element model. A cycle jump means that the com- 
putation is done for a certain set of loading cycles at chosen intervals, and 
that the effect on the stiffness degradation of these loading cycles is extrap- 
olated over the corresponding intervals in an appropriate manner. The cycle 
jump strategy adopted here is based on the one presented in [36]. After a 
certain number of cycles N l} the damage variable d{ at an integration point 
J is computed using the quasi-static constitutive equations. The predicted 
evolution of the damage variable with the number of cycles, J^, is computed 
using Equation (19). The damage variable at an integration point J after ANi 
cycles is: 


-\f 

h+ANi 


= d + 


ddJ 

— AA^- 
ON 4 


(30) 


To determine the number of cycles ANi that can be skipped with a controlled 
level of accuracy, the following equation is used: 


ANi 


Ad max 



r 

max < 
j 

iW 


(31) 


where Ad max is a pre-established value. The smaller the choice of Ad max the 
higher the accuracy of the analysis. 
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5 Results and discussion 


The present model is implemented as a user-written finite element in ABAQUS® 
[37] by adding the fatigue damage model to the constitutive behavior of a co- 
hesive element previously developed in Refs. [13,14]. 

Several single-element tests were performed to verify the response of the fa- 
tigue damage model. Then, simulations of Mode I, Mode II and mixed-mode 
delamination tests were conducted to demonstrate that when the constitutive 
damage model is used in a structural analysis, the analysis can reproduce the 
response of the test specimens without the use of any model-specific adjust- 
ment parameters. 


5.1 One element tests 


The finite element model shown in Figure 7 is composed of two 4- node plane 
strain elements connected by a 4-node cohesive element representing the in- 
terface. 


◄ — 

Applied 

Displacement 




Fig. 7. Undeformed mesh with the boundary conditions and deformed mesh of one 
cohesive element tests. 


The material properties shown in Table 2 correspond to a T300/977-2 carbon- 
fiber reinforced epoxy laminate. The Paris Law coefficients used in the simu- 
lation were C\ = 0.0616 mm 2 /cycle and mi = 5.4. The threshold for fatigue 
crack propagation was assumed to be zero. 


The load was applied in two steps. The first step was a quasi-static step where 
the displacement jump is incremented to 20 times the onset displacement. The 
second step accounts for fatigue damage resulting from a maximum applied 
displacement of 20 times the onset displacement and a load ratio R = 0. 
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Table 2 

Properties used in the models with only one cohesive element. 


En 

E' 2'2 — E 33 

G 12 — G 1 3 

G 23 

Ul2 = U13 

U23 

150.0 GPa 

11.0 GPa 

6.0 GPa 

3.7 GPa 

0.25 

0.45 

Gic 

Giic 

T° 

r 3 

T° — T° 
'2 _ T 1 

K 


0.268kJ/m 2 

0.632 kJ/m 2 

45 MPa 

45 MPa 

10 6 N/mm 3 




Fig. 8. Evolution of the interface traction in the constitutive equation for a displace- 
ment jump controlled high-cycle fatigue test. 


The evolution of the interface traction in the constitutive equation for a high- 
cycle fatigue test under displacement control is shown in Figure 8. It can be 
observed that fatigue damage causes a reduction of the stiffness, the interfacial 
traction, and the interfacial strength. The evolution of the interface traction 
and strength with the number of cycles is shown in Figure 9. The shape of the 
obtained curves is similar to the widely-used S-N curves used in the design for 
fatigue strength. 


5.2 Simulation of a DCB specimen under fatigue loading 


Simulations of a double-cantilever beam (DCB) specimen were conducted to 
simulate the crack growth rate under Mode I loading for different ranges of the 
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10° 10 1 io 2 io 3 io 4 10 5 10° 

Cycles 


Fig. 9. Evolution of the interface traction and the maximum interface strength as 
a function of the number of cycles for a displacement jump controlled high-cycle 
fatigue test. 


energy release rates. Experimental data on fatigue-driven delamination growth 
reported by Asp et al. [38] was selected for the validation of the numerical 
model. The specimen was fabricated with HTA/6376C carbon/epoxy prepreg 
produced by Hexcel. The layup consisted in [0i 2 //(±5/0 4 )s], where the sign 
/ / refers to the plane of the artificial delamination. The specimen was 150- 
mm-long, 20.0 nun-wide, with two 1.55-mm-thick arms, and an initial crack 
length of 35nnn. A description of the experimental procedure is reported by 
Asp et al. [38]. The material properties are shown in Table 3 [11,38,39]. 


Table 3 

Material properties for HTA/6376C carbon/epoxy [11,38,39]. 


An(GPa) 

E 22 — A 33 (GPa) 

G12 = Gi 3 (GPa) 

G 2 3 (GPa) u \2 — v\ 3 u 2 3 

120 

10.5 

5.25 

3.48 0.30 0.51 

Gi c (kJ/m 2 ) 

Gii c (kJ/m 2 ) 

rij 1 = rf (MPa) 

T 3 (MPa) 

0.260 

1.002 

30 

30 
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In the finite element model, the specimen’s arms are loaded with opposing 
moments (Figure 10) to obtain a Mode I energy release rate that is independent 
independent of crack length and, consequently, to achieve a constant fatigue 
crack growth rate. The energy release rate is related to the applied moment 
as [11]: 


G\ = 


M 2 

bEI 


(32) 


where b is the specimen width, E is the longitudinal flexural Young’s modulus 
and / is the second moment of area of the specimen’s arm. 



The finite element model is composed of 4-node plane strain elements for 
the arms, which are connected by 4-node cohesive elements representing the 
interface. Two elements are used through the thickness, h, of each arm. The 
length of the element is 0.05nnn (see Figure 11). The Paris Law parameters of 
Equation (1) were obtained from a linear regression of the experimental data 
[38]: C\ = 0.0616 mm 2 /cycle and up = 5.4. The energy release rate threshold 
is 0.060 kJ/m 2 [38]. The material properties are shown in Table 4. 


Initial crack 
(35mm) 




Fig. 11. Detail of the FEM model of the DCB specimen. Two applied load P with 
opposite direction were applied to each arm. The applied moment is equal to the 
product between the applied load P and the thickness of the arm. 


The load is applied in two steps: the first analysis loading step is quasi-static 
and it ends at the maximum applied load. It is assumed that no fatigue damage 
accumulation occurs during this step. Next, a second loading step is applied 
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in which the maximum load is held constant, during this cycle, the analysis 
pseudo-time increment is assumed to be proportional to the number of loading 
cycles so that the fatigue damage model accounts for the accumulation of cyclic 
damage. The maximum variation in the damage variable Ad max allowed in a 
cycle jump is set to 0.001. 

Table 4 

Fatigue material properties for HTA/6376C carbon/epoxy obtained from references 
[11,39] and using Equations (27) to (29). 


Ci (mm/cycle) 

C\\ (mm/cycle) 

CWo (mm/cycle) 

m m (mm/cycle) 

0.0616 

2.99 

4.23 

458087 

mi 

mu 

m 50% 

Cm 

5.4 

4.5 

6.41 

4.94 

Gith (kJ/m 2 ) 

Gmh (kJ/m 2 ) 

^U30%th (kJ/m 2 ) 

V 

0.060 

0.100 

0.066 

2.73 


The results obtained from the simulations and the experimental data are 
shown in Figure 12. It can be observed that the constitutive model accounts 
for all three regions of fatigue crack growth. In region II, where crack growth 
rates follow the Paris Law, it is observed that a good agreement between the 
predictions and the experimental data is obtained. In region I there is negli- 
gible crack growth rate for small values of the normalized energy release rate 
and the numerical data follows the trend of the experimental data. A signifi- 
cant difference between the numerical and the experimental data is observed 
in region III. One of the reasons for this difference is that the crack growth 
rates present in region III are very high and, therefore, a low-cycle instead 
of a high-cycle fatigue model is more appropriate for this region. However, in 
spite of this difference, the model can also predict Region III crack growth 
rate, where the Paris Law equation is not valid. 
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Fig. 12. Comparison of the experimental data with the crack growth rate obtained 
from the numerical simulation for a Mode I DCB test. 


5.2.1 Sensitivity of the rack propagation rates to the load ratio 

Several DCB tests were conducted to verify the sensitivity of the model to 
the load ratio. The results obtained from the simulations are shown in Figure 
13 where it can be observed that higher load ratios decrease the crack growth 
rate. 

The sensitivity of the constitutive model to the load ratio is an asset of the 
model. The sensitivity of the propagation rate to the load ratio derives directly 
from the quasi-static model rather than from a fatigue model defined as a 
function of the load ratio, as has been done in previous investigations [12]. 

It can be observed from Figure 13 that the same energy release rate threshold 
G t h is predicted for all load ratios. This result is a consequence of the current 
formulation of the model: the influence of the load ratio on the energy release 
rate threshold is not taken into account. If G t h is constant, then the energy 
release rate range threshold A G t h must vary with the load ratio, which is 
a trend not reflected in experimental results. To verify this dependence, the 
predicted crack growth rates are represented in Figure 14 as a function of the 
energy release rate range AG instead of the G max . It can be observed that the 
predicted crack growth rates for the Region II of the crack growth rate are 
almost independent of the load ratio, which is in agreement with experimental 
results [40]. However, it can also be observed that different energy release 
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Fig. 13. Sensitivity of the model to the load ratio for a Mode I DCB test. 

rate range thresholds are predicted for different load ratios. This effect is a 
limitation of the present model that will be addressed in future work by the 
use of the range of the energy release rate range threshold in Equation (21). 



Normalized Energy Release Rate Range, A G/G h 


Fig. 14. Sensitivity of the model to the load ratio for a Mode I DCB test. 
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5.3 Simulation of a fENF test 


Several tests were conducted to simulate the crack growth rate under Mode II 
loading for different ranges of the energy release rate. Experimental data on 
fatigue driven delamination growth reported in [38] was selected for compari- 
son. The dimensions and the material of the specimen are the same used for 
the DCB specimen described in the previous section. 

For pure Mode II, the specimen was loaded using the four point End Notched 
Flexure (4ENF) test shown in Figure 15. The energy release rate is related to 
the applied moment, f , as [11]: 


3(f) 2 
4 bEI 


(33) 



The hnite element model used was similar to that used in the simulation of the 
Mode I test (see Figure 16). The Paris Law coefficients of Equation (1) were 
obtained from a linear regression of the experimental data presented in Ref. 
[38]: C\\ = 2.99 mm 2 /cycle and mu = 4.5. The energy release rate threshold 
is 0.100 kJ/m 2 [38]. The fatigue properties are summarized in Table 4. 

The load is applied in two steps, as described in the previous section. 

The results obtained from the simulations and the experimental data [38] 
are shown in Figure 17. The predicted crack growth rates for small values of 

gmax 

-gh— are slightly higher compared to the experimental data. This difference 
can be attributed to friction effects that are not considered in the current 
implementation. Moreover, it should be mentioned that the model estimates 
the size of the cohesive zone using Equation (20), which may be accurate for 
Mode I [33] but not necessarily for Mode II. Further investigations on the 
estimation of the cohesive zone length under Mode II should be conducted. 
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c 


p 



p 



Fig. 16. Detail of the FEM model of the fENF specimen. 



Normalized Energy Release Rate, Gff*/G Jh 


Fig. 17. Comparison of the experimental data with the crack growth rate obtained 
from the numerical simulation for a Mode II fENF test. 
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5-4 Simulation of mixed-mode loading 


Several tests were conducted to simulate the crack growth rate under mixed- 
mode loading with Gi = Gu for different energy release rates. Experimental 
data on fatigue driven delamination growth reported in [38] was selected for 
comparison. The dimensions and the material of the specimen are the same 
used for the DCB specimen described above. 

For mixed-mode loading, the specimen was loaded with two moments as it 
is shown in Figure 18. The ratio between the two applied moments, p, for a 
mode-ratio of 50% is: 


P 



(34) 


m Ci 
pM 0 




Fig. 18. Loading pattern for mixed-mode specimen. 


The energy release rate is related to the applied moment, M , as [11]: 


G\ — Gu — 


M- 


i(i + f) bEI 


(35) 


The finite element model used was similar to that used in the simulation of 
the Mode 1 test (see Figure 19). 


The mixed-mode parameters C m , m m , and G t h are computed at each inte- 
gration point using Equations (27), (28) and (29) to account for any changes 
in the mode-ratio. The fatigue material properties used in the simulation are 
summarized in Table 4. The load is applied in two loading steps, as described 
in previous sections. 

The results obtained from the simulations and the experimental data [38] are 
shown in Figure 20. As in the case of pure Mode II, the predicted data for 
small values of are slightly higher than the experimental data. 
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Fig. 19. Detail of the FEM model of the specimen mixed-mode loaded. 



Fig. 20. Comparison of the experimental data with the crack growth rate obtained 
from the numerical simulation for a mixed-mode test with Gj = Gn- 
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6 Conclusions 


A damage model suitable for both quasi-static and high-cycle fatigue delami- 
nation propagation was developed. The evolution of the damage variable was 
derived by linking Damage Mechanics and Fracture Mechanics, thus establish- 
ing a relation between damage evolution and crack growth rates. The damage 
evolution laws for cyclic fatigue were combined with the law of damage evolu- 
tion for quasi-static loads within a cohesive element previously developed by 
the authors. 

The model was validated using single-element numerical tests, as well as by 
simulating the propagation rates of Mode I, Mode II and mixed-mode tests. 
The model was able to reproduce the Paris Law growth rate without the need 
of any additional adjustment parameters. Moreover, the model accounts for the 
energy release rate thresholds preventing crack growth for smaller values of the 
energy release rate. Unlike other approaches proposed in the literature, where 
the dependence on the load ratio, R , is introduced through the definition of in- 
dependent Paris Law parameters, the effects of the load ratio on the analysis 
results is inherent to the formulation. 

The analysis of the results indicate that the model is more accurate when 
Mode I loading predominates. This effect can be justified by two factors: i) 
friction between the crack faces is not taken into account in the model, and ii) 
the equation used to estimate the length of the cohesive zone was developed 
for for pure Mode I loading. Further investigations on the estimation of the 
cohesive zone length under Mode II and mixed-mode should be conducted. 

In summary, the model is able to predict the crack growth rates in all regimes 
of propagation and the results compare favorably with the experimental data, 
including the negligible crack growth rates for small values of the normalized 
energy release rate and the sensitivity to the mode and load ratio. 
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